Energy transport in jammed sphere packings 
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We calculate the normal modes of vibration in jammed sphere packings to obtain the energy 
diffusivity, a spectral measure of transport. At the boson peak frequency, we find an Ioffe-Regel 
crossover from a diffusivity that drops rapidly with frequency to one that is nearly frequency- 
independent. This crossover frequency shifts to zero as the system is decompressed towards the 
jamming transition, providing unambiguous evidence of a regime in frequency of nearly constant 
diffusivity. Such a regime, postulated to exist in glasses to explain the temperature dependence of 
the thermal conductivity, therefore appears to arise from properties of the jamming transition. 

PACS numbers: 45.70.-n, 61.43.Fs, 65.60.+a, 83.80.Fg 



Zero-temperature soft-sphere models inspired by foams 
and granular media have given insight not only into the 
geometry of hard-sphere packings [J, d, 0| but also into 
the physics of the low-energy excitations in glasses [U, Q . 
In particular, a model system of frictionless spheres inter- 
acting with finite-ranged repulsions exhibits a jamming 
transition (Point J) at a density corresponding to ran- 
dom close-packing of hard spheres [iL At the transi- 
tion, the coordination number jumps l], [|| from zero to 
the minimum value required for mechanical stability, the 
"isostatic" value Q, and there is a plateau in the den- 
sity of vibrational states that extends down to zero fre- 
quency [|(|. Upon compression, this plateau persists but 
only above a characteristic frequency, uj*, that increases 
with density. The modes in the plateau region have been 
shown to arise from zero-frequency vibrational modes 
at the isostatic transition 0]. These anomalous modes 
are in excess of the Debye prediction and are directly 
connected d, Q to excess vibrational modes in glasses, 
known as the "boson peak" (Til. Il2l Il3j . However, it is 
not clear how these modes contribute to heat conduction. 

In this paper, we investigate thermal transport as a 
function of compression in jammed sphere packings. At 
the jamming threshold, we find that all delocalized modes 
transport heat with a low diffusivity nearly independent 
of frequency, in contrast to ordinary solids in which sound 
modes transport heat ballistically with a diverging diffu- 
sivity in the long- wavelength limit. The behavior at Point 
J is reminiscent of many amorphous solids, which unlike 
crystals display a thermal conductivity that rises mono- 
tonically with temperature T |14| . This property has 
been posited to arise from a frequency regime of small, 
constant diffusivity [1(1 HE Gil- We show that this regime 
can originate from the vibrational spectrum at Point J. 
Upon compression, the low-diffusivity modes persist, but 
only above a crossover frequency corresponding to the 
frequency of the boson peak, lo* [l7l [Lsj ; below ui* , the 
spectrum is dominated by transverse plane waves. 

Our model [l], H[ is a 50/50 mixture of frictionless 
spherical particles with a diameter ratio of 1.4. Parti- 
cles i and j interact in three dimensions via a one-sided 



harmonic potential: U (r^ 
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distance between their centers, r 
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— fij/cij) 2 when the 
is less than the sum 
of their radii, cr^ and zero otherwise. Jammed pack- 
ings at temperature T — are obtained by conjugate- 
gradient energy minimization. We study systems of 
250 < N < 10, 000 particles with periodic boundary con- 
ditions in all directions. The packing fraction at the on- 
set of jamming, <p c , is characterized by the onset of a 
nonzero pressure. We determine <j) c and obtain T = 
configurations at controlled A0 = cf> — 4> c as in Ref. [1] . 
For each configuration, we diagonalize its Hessian matrix, 
whose m th eigenvalue is the squared frequency, w^, of the 
orthonormal eigenmode described by the displacement 
e m (j) of each particle j. We used the package ARPACK 
to handle large systems [l9j] . The particle mass, M, in- 
teraction energy, e, and diameter of the smaller particle, 
<7, are set to unity. The frequency is in units of \J 'e/Ma 2 . 

We also study an "unstressed" model in which we use 
energy-minimized configurations obtained from the pre- 
vious model, and replace the interaction potential U(rij) 
between each pair of overlapping particles with an un- 
stretched spring with the same stiffness, U"(rij). Be- 
cause all springs are unstretched, there are no forces 
between particles in their equilibrium positions so that 
stable configurations for the stressed system are also sta- 
ble in the unstressed one. This corresponds to dropping 
terms depending on U' in the Hessian @, [1] . 

For a strongly scattering system, a diffusive description 
of energy transport can be more useful than one in terms 
of ballistic propagation with a very short mean free path. 
Therefore instead of calculating the thermal conductiv- 
ity, k(T), directly using molecular dynamics HE [SI, we 
calculate the thermal diffusivity, d(ui m ) for vibrational 
mode 771 [IE [Uj]. k(T) can be expressed in terms of 
d(uj m ) and the heat capacity C{ui m ) (TBI. [HI |22| : 

1 1 f°° 

k(T) - - ]T C{uo m )d{uj m ) = - du;D(Lj)C(Lj)d(Lj), 

m 

(i) 

where the sum runs over all vibrational modes m, V is 
the volume of the system, and D(u>) = Yl m H^m ~ is 
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the density of vibrational states. Thus, in Eq. [TJ C(ui) = 
keifihjjfe^ /{e phuJ - l) 2 (where (3 = l/k B T and k B is 
the Boltzmann constant) depends on T. It characterizes 
the heat carried at frequency u>, while d(u>), which has 
units of <jy/ e/M, is a T-independent scattering function. 

The physical meaning of the diffusivity is best illus- 
trated operationally. Consider a wavepacket narrowly 
peaked at to and localized at r at time t = 0. Over time, 
the wavepacket spreads out and can be characterized by 
a time-independent diffusivity given by the square of the 
width of the wavepacket at time t, divided by t [1 51 ] - In a 
weakly-scattering system, d(uj) = cl(u))/3, where c is the 
sound speed and t{ui) the phonon mean-free path. 

We follow Allen and Feldman [l6| to calculate d{u>) 
in terms of the normal modes of a given configuration, 
using the Kubo-Greenwood formula for the thermal con- 
ductivity as the response to a temperature gradient that 
couples different modes. We use [la ] 
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where the vector heat-flux matrix elements are 



|S(c,c')| 2 = 



Z)mn \£"mn\ 2 8(u - U m )8(bj' - LJ n ) 

D(uj)D(u)') 



(3) 



where m and n index the vibrational modes. 

For a finite system the modes are discrete. We calcu- 
late the matrix elements S mn from the Hessian H % K and 
its m th normalized eigenvector e m (i;a) [161 ] via 

S mn = 2J (*t - rj)e m (i; a) e n (j; (3), (4) 

where {i, j} and {a, /?} label particles and their Cartesian 
coordinates respectively. 

In a finite system, the delta function in Eq. [2] must 
be replaced by a representation with nonzero width, 77. 
We use [HI g(uj m - ui n , 77) = rj/[n((u m - u; n ) 2 + if )] and 
define 
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We set 1] = jSu), where 7 > 1 and Slo is the average 
spacing between successive modes. The desired d(u)) is 
then d(u>, r], N) in the double limit N — > 00, 77 — > + . 

We now use Eqs. [2-5] to calculate the diffusivity. Our 
goal is to extract d(tu) for an infinite system so we must 
confront finite-size effects. We will show that |S| 2 has a 
particularly simple form at the jamming threshold, en- 
abling us to determine the N — > 00 behavior. 

Figure HJa) shows the heat-flux matrix elements 
|£(u>, u')\ 2 defined in Eq.[3]for packings at <p~<p c = 10~ 6 
for different values of oj versus ui'. Fig.QJb) shows that all 
the curves can be collapsed for different system sizes, JV, 
and frequencies except at high lu' where the modes be- 
come localized HH HH ■ The inset to Fig.[TJb) shows 



FIG. 1: (Color Online) Diffusivity just above the jamming 
transition at A(j> — 10 -6 . (a) Heat-flux matrix elements 
|E(oj,tj')| 2 plotted versus u)' at N = 2000 for lu = 0.012 
(solid), 0.035 (long dashed), 0.08 (short dashed), 0.10 (dot- 
ted), 0.22 (dot-dashed), and 0.68 (dot-dash-dashed), (b) Scal- 
ing plot showing collapse of \T,(lu,lo ')| 2 at TV = 2000 (black 
solid), 1000 (red dashed), and 500 (blue dotted) with scale 
factors s 2 and w. The green dot-dashed line indicates a power- 
law slope of 2. Inset: Scale factors s 2 (red symbols) and w 
(black symbols) versus lu. We find s 2 oc lu 2 (red dashed line) 
and w <x lu (black dotted line) except at high lu. (c) Plot of 
d(Lu,rj = 0.002,/V = 2000) defined in Eq. © (dashed). Solid 
line: predicted d(tu) for TV — > 00. 



that the scale factors for the collapse satisfy s 2 — lu 2 and 
w = u>, respectively, except at high frequencies where lo- 
calization sets in. Note that the scaling collapse demon- 
strates that the only noticeable system-size dependence 
is a prefactor of 1//V- Since for large N, the density of 
77)states scales as N Q, Eq. [^therefore yields a well-defined 
diffusivity in the N — > 00 limit (solid curve in Fig. [He)). 

Note that the collapse in Fig. QJb) implies that 
E(cj,cj)| 2 oc Qj 2 /N a.t low frequencies. This scaling arises 
when overlap with nearby modes, described by Eq. ((¥]), 
is small and independent of frequency and when modes 
are spatially uncorrelated [2(| . Thus Eq. [2] implies that 
d(u) oc D(lo) at low u>. At Point J, because the density 
of states is nearly constant down to ui — 0], the diffu- 
sivity is nearly constant as well. (The small slope of d(ui) 
with frequency is due to the slight frequency dependence 
of D(ui).) These results show that although the low- 
frequency modes are extended at the jamming thresh- 
old [24|], they do not behave like plane waves and the 
usual divergence of diffusivity is completely suppressed. 

Over most of the frequency range, this N — > 00 predic- 
tion agrees very well with the dashed curve in Fig. [IJc), 
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FIG. 2: (Color Online) Diffusivity versus compression. 
d(to,r) = 0.002, N = 2000) for the (a) unstressed and (b) 
stressed systems at A<f> — 10 -6 (solid black), 0.01 (red 
dashed), 0.05 (green dotted), 0.1 (blue dot-dashed), and 0.3 
(purple dot-dot-dashed). In (a) the cyan dashed line indicates 
a power- law of uj~ 4 and the closed symbols indicate the degen- 
erate sets of discrete plane- wave modes in our finite system. 
Inset: The ratios Ud/io* (open circles) and tOd/tom (open tri- 
angles) versus A(j>. (c) d(io,r] — 0.004, N) for the stressed 
system at A(j> = 0.5 for N = 10,000 (black solid), 2000 (red 
dashed), and 500 (blue dot-dashed). 



which shows d(to, r] — 0.002, N — 2000) for a finite sys- 
tem. However, at low frequency d(ui,r),N) exhibits an 
upturn; this upturn is a finite-size artifact that can be 
shown from Eqs. [2-5] to scale as cu~ 3 with a prefactor 
that vanishes as N — > oo, 77 — ^ [26]. 

We now study diffusivity as the system is compressed 
beyond the jamming threshold. We begin with the un- 
stressed model where the data are particularly simple to 
interpret. In this case, the system is always held at zero 
pressure, so that increasing A<ft = <f> — <j) c corresponds 
to increasing only the average coordination number of 
the network of interacting particles. At all compressions, 
Fig. [Ha) shows that d(to, rj = 0.002, N = 2000) vanishes 
at high to for localized modes. At low frequencies, there is 
the upturn due to finite-size effects discussed above; this 
upturn lies at frequencies below those shown in Fig. [21 At 
intermediate frequencies, the diffusivity decreases rapidly 
with increasing lo until it reaches a small constant value 
do for to > u>d- Below ui^, the modes are discrete due 
to the finite system size, as indicated by the discrete 
points in Fig. [2ja). The mode frequencies correspond 
to lu = cxkn, where ct is the transverse sound speed and 
k n are the lowest allowed wavevectors. The calculated 
diffusivity below cod decreases sharply with increasing lu 



as expected for scattering of plane waves [221 ]. Thus, u>d 
marks the crossover from transverse plane waves to a 
small, nearly constant diffusivity. 

The frequency iOd can be understood as the Ioffe-Regel 
crossover from weak to strong scattering of transverse 
modes [27]. For lu < LUd, the transverse plane-waves obey 
lu = Cxk, where k is the wavevector. As lu approaches LUd, 
we have shown d(u>) = ct£/3 —> do, where £ is the phonon 
scattering mean free path. The Ioffe-Regel criterion, kl w 
1, predicts a crossover frequency near lujr = c T /3do. For 
our harmonic system, the transverse speed ct oc Ac/) ' 25 
Q so that loir oc Acj) - 5 . The inset to Fig. H indeed 
shows for transverse modes, LOd/iOiR ~ 2 over our range 
of compression, consistent with the Ioffe-Regel criterion. 

Fig. [^a) (inset) also shows the ratio of LUd to the bo- 
son peak frequency, to*, defined as the onset frequency for 
the plateau in the density of states, which was previously 
shown to scale as A</> 5 0, The ratio LUd/co* is con- 
stant over a wide range of A(j>. Studies of silica [H, [29[ 
and several disordered models [lj], OH observe that the 
boson peak frequency and Ioffe-Regel crossover frequency 
agree within a factor of order unity. In our system, this 
relation is unambiguous because both frequencies shift 
together as Acf> is varied. Our results also indicate that 
at Point J, where lu* — 0, the diffusivity should remain 
nearly flat down to zero frequency, as shown indepen- 
dently in Fig.QJc). Therefore, the modes above LUd (i.e., 
those with constant diffusivity) can be identified as the 
anomalous modes that derive from soft modes at the iso- 
static point [1, [§] . 

While the unstressed system may be appropriate for 
systems where the coordination at threshold exceeds the 
isostatic value (such as frictional systems [3(|), we are 
also interested in systems where inter-particle forces in- 
crease with Acf>. These forces lower both the sound speed 
ct which controls the Ioffe-Regel crossover frequency [lj 
and the frequencies of modes in the plateau p|. Finite- 
size effects, which cut off plane waves at low lu, are there- 
fore more obstructive in stressed systems. 

Fig. [Hb) shows that there is no discernible change 
in the diffusivity of the stressed system over the range 
10~ 6 < A(j) < 10~ 2 . Above A0 « 10~ 2 , structure de- 
velops at intermediate frequencies. Each peak can be 
identified with one of the first few allowed wavevectors 
for longitudinal or transverse modes [2(|. Fig.[^c) shows 
that as N increases, the plane-wave peaks shift to lower 
frequencies and grow closer together, as expected, so that 
for high enough N, peaks in a given frequency range will 
merge into a smooth curve as peak widths exceed their 
spacing. We therefore conclude that this structure will 
disappear in the infinite-size limit when the only observ- 
able plane-wave peaks will be shifted to zero frequency. 

Fig- (He) also shows that the peak heights increase with 
N and decreasing lu at low frequency. This suggests that 
the diffusivity rises smoothly with decreasing lu at fre- 
quencies below some LUd in the thermodynamic limit, and 
has a constant value do above 10 d- In the unstressed case, 
we found plane-wave behavior below LUd ~ loir — c T /3do. 
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In the stressed case, we speculate that u/4 might be sim- 
ilarly defined and should also increase with A(j) - 5 since 
ct oc A<fi°- 25 for the stressed as well as the unstressed 
case 

At the jamming threshold, we can calculate the ther- 
mal conductivity, n{T), from Eq. [T] and obtain a finite 
answer even within the harmonic approximation because 
the diffusivity does not diverge at low oj. By setting d(uj) 
and D(oj) to be constant up to the localization thresh- 
old, we find k oc T up to the Debye temperature. At 
high temperatures, k saturates when all the modes are 
excited. 

One of the most striking differences between heat con- 
duction in ordered and disordered structures is that 
the thermal conductivity n of crystalline materials first 
rises with increasing T but eventually drops due to 
phonon-phonon scattering [14| , while n for glasses in- 
creases monotonically in T. This perplexing property 
of glasses has been explained heuristically by assuming 
that phonons are scattered so strongly by disorder that 
transport becomes diffusive, with a frequency regime of 
small, constant thermal diffusivity. In that case, the 
thermal conductivity simply increases with the heat ca- 
pacity according to Eq. Q] [l(| HE [IB]. In °ur finite- 
sized unstressed systems, we see clear evidence for such 
a regime of nearly constant diffusivity, and find that its 
frequency onset, given by the Ioffe-Regel crossover for 



transverse phonons, increases with compression. For the 
compressed system with stress, our results are much less 
clear in the low-w regime because of finite-size effects. 
However, above some intermediate frequency, we again 
find a constant diffusivity and constant density of states, 
which lead to the rise in the thermal conductivity as in 
the unstressed case. 

In earlier work, we showed that the vibrational spectra 
of model systems such as the Lennard- Jones glass could 
be understood in terms of jammed sphere packings [j|. 
We now find that these packings capture some of the 
crucial physics invoked to explain the temperature de- 
pendence of the thermal conductivity-a crossover from 
low-cj transverse phonons to excess vibrational modes of 
nearly constant diffusivity [13] • The physical origin of 
the diffusive modes lies in the behavior of packings at 
the jamming threshold, where the Ioffe-Regel crossover 
frequency vanishes. Upon compression, the flat diffusiv- 
ity shifts to higher frequencies but does not disappear. 
Thus, compressed sphere packings are a useful starting 
point for understanding energy transport in glasses. 
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